Slow Switching in Globally Coupled Oscillators: 
Robustness and Occurrence through Delayed Coupling 
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The phenomenon of slow switching in populations of globally coupled oscillators is discussed. This 
characteristic collective dynamics, which was first discovered in a particular class of the phase 
oscillator model, is a result of the formation of a heteroclinic loop connecting a pair of clustered states 
of the population. We argue that the same behavior can arise in a wider class of oscillator models 
with the amplitude degree of freedom. We also argue how such heteroclinic loops arise inevitably 
and persist robustly in a homogeneous population of globally coupled oscillators. Although the 
heteroclinic loop might seem to arise only exceptionally, we find that it appears rather easily by 
introducing the time-delay in the population which would otherwise exhibit perfect phase synchrony. 
We argue that the appearance of the heteroclinic loop induced by the delayed coupling is then 
characterized by transcritical and saddle-node bifurcations. Slow switching arises when the system 
with a heteroclinic loop is weakly perturbed. This will be demonstrated with a vector model by 
applying weak noises. Other types of weak symmetry-breaking perturbations can also cause slow 
switching. 
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I. INTRODUCTION 

Coupled limit-cycle oscillators appear in various con- 
texts in physics jjj-U, chemistry |||| and biology [@ |). 
Various types of collective behavior which arise when 
they form large assemblies have been studied extensively 
over the last few decades. Among the possible types of 
behavior, we will particularly be concerned with cluster- 
ing and slow switching which was first studied by Hansel 
et al. @ in a homogeneous population of globally cou- 
pled phase oscillators. Assuming a simple form for the 
coupling function, they showed numerically that after a 
long transient the system approaches a two-cluster state, 
i.e. the whole population splits into rigidly rotating two 
subpopulations each in perfect phase synchrony. How- 
ever, the stability analysis of this two-cluster state re- 
vealed that it is linearly unstable, corresponding to a 
saddle point if seen in a co-rotating frame of reference. 
The seeming contradiction here was interpreted in terms 
of the formation of a heteroclinic loop connecting this 
two-cluster state and another two-cluster state which is 
obtained simply by a constant phase shift of the former. 
In fact, when this heteroclinic loop is attracting, the tra- 
jectory comes to stay longer and longer near these saddle 
points, so that the numerical round-off error finally forces 
the system to stay at one of the saddles forever. This in- 
terpretation was supported by the fact that when small 
external noise is included the system is no longer fixed 
at a saddle point but starts to repeat slow switchings be- 
tween the pair of saddles (see Fig. [j]). Although these 
findings are so important, explanation of the reason is 
still needed why the heteroclinic loop arises inevitably 
and persists robustly against our common belief in its 
structural instability. 

In the next two sections, we restrict our consideration 
to the phase model. In Sec. ||, we discuss in some detail 
the stability condition of the two-cluster state. Existing 
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FIG. 1. Slow switching exhibited by the model in [10]. The 
figure displays the time-evolution of the number density of the 
oscillators as a function of the phase. The whole population, 
which was initially almost uniform, splits into two subpopu- 
lations each almost converging to a point cluster. After some 
time, however, this seeming convergence turns out unstable, 
and is followed by the period of their scattering, but again 
this is followed by the period of the convergence, and so forth. 
This form of alternation between the two characteristic peri- 
ods of the convergence and dispersion of the clusters is called 
slow switching. The phase-advanced/retarded cluster at the 
end of one cycle becomes phase-retarded/advanced at the end 
of the next cycle. 

numerical results obtained by a particular model sug- 
gest the appearance of heteroclinic loops. Thus, in Sec. 
Ill , we argue the mechanism by which heteroclinic loops 
are necessarily formed. Specifically, a sufficient condition 
will be given for the existence of a heteroclinic loop, and 
how this condition is satisfied in the phase model will 
be discussed. In Sec. IV, we introduce a specific vec- 



tor oscillator model for globally coupled oscillators and 
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exhibit numerically that heteroclinic loops are formed in 
our vector model. We show there that the phase-coupling 
function, derived numerically from the vector model by 
the method of the phase reduction, satisfies the above- 
mentioned condition leading to the formation of hetero- 
clinic loops. Our result gives the first result showing the 
heteroclinic loop in the vector model of coupled oscilla- 
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tors. In Sec. |V|, we generalize the argument in Sec 
to the vector model. 

The formation of heteroclinic loops in globally coupled 
oscillators may seem to be a pathological phenomenon 
which can happen only exceptionally. However, the time- 
delay in coupling can easily cause a bifurcation from per- 
fect synchrony to the formation of a heteroclinic loop, 
and this will be discussed in Sec. VI. The method of the 



phase reduction provides a clear understanding of why 
this is actually possible. Slow switching becomes per- 
sistent when the oscillators are subject t o we ak external 
noise, which will be discussed in Sec. VII by using a 
vector model. We will also show there that the same 
phenomenon can also be caused by other types of ran- 
domness. 



II. HETEROCLINIC LOOP IN THE PHASE 
MODEL 

Populations of weakly coupled limit cycle oscillators 
can be described by the phase model p[. Suppose that 
the oscillators are identical, each interacting with all 
the others with equal strength. Then the corresponding 
phase model is expressed as 



dt 



K 

N 



N 



(i) 



where "4>i(t) (0 < ipi < 2ir) is the phase of the i-th oscil- 
lator (i — 1, • • • , N), u> and K are positive constants, and 
T(x) is the coupling function with 27r-periodicity. 
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FIG. 2. Long transient of the order parameter after which 
the whole population converges to a two-cluster state. 
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FIG. 3. Condition for the existence of two-cluster 
states. A takes three values for a given p within the range 

Pmin < p < 1 - Pmin- 

Hansel et al. [[l0| analyzed the case of a particular form 
of the coupling function 



T(x) = - sin(x + 1.25) + 0.25 sin(2x) 



(2) 



They showed by numerical simulations that the oscilla- 
tors with random initial distribution are assembled to 
form two subgroups each in perfect phase synchrony but 
with a constant mutual phase difference. The collective 
behavior of the system can conveniently be described in 
terms of the order parameter defined by 



O(t) 



1 

N 



N 



(3) 



Its value is 1 for perfect synchrony and for perfect in- 
coherence. A time trace of the order parameter for the 
above model is displayed in Fig. ||. The oscillators be- 
longing to the respective groups are identical in phase, 
and this pair of point clusters rotate rigidly at a constant 
angular frequency. The mutual phase difference is de- 
noted by A. We call hereafter the phase-advanced and re- 
tarded clusters the A-cluster and B-cluster, respectively. 
Let the fraction of the oscillators belonging to the A- 
cluster be p. Such a two-cluster state may thus be speci- 
fied by (p, A), where A is within the region — tt < A < 7r 
by convention. This set of values may generally differ for 
different initial conditions. 

Existence and stability of the two-cluster states can be 
analyzed as follows. Consider a two-cluster state with 
phases -0 a and ip-B. Eq. (pi) then becomes 



dt^ 



w + K{pT(0) + (l-p)r(x)}, 



(4) 



-Mt) = u + K{(1 - p)T(0) + pT(-x)}, (5) 

where x denotes the phase difference, i.e. x = -0a — ^b- 
Subtracting Eq. (H) from Eq. (0), we obtain 
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FIG. 4. Eigenvalues about two-cluster states as a function 
of A. All two-cluster states are unstable here. Note that the 
eigenvalues A2 and A3 are negative for the states (p, A') and 
(1-P,A"). 



-x{t) = K{(2p - l)r(O) + (1 - p)T(x) - pT(-x)}. (6) 
Since x is constant in the two-cluster state, we have 



P (A) 



r(o) - r(A) 



2r(o)-r(A) -r(-A)' 



(7) 



(p, A) exists with p satisfying < p < 1. Substituting 
Eq. (|) into Eq. we obtain the condition for the 

existence of (p, A), and this is displayed graphically in 
Fig. H A takes three values for a given p within the 
range Pmln < p < 1 — Pmin j where Pmin is defined by the 
minimum value of p in the range < A < n. These three 
states are denoted by (p, A'), (p, - A") = (1 -p, A") and 
(p, A'"), where A' and A" are understood to be positive 
and j A'" I to be larger than A' and A". 

The eigenvalues of the stability matrix are given by 

A = 0, (8) 
Ai = K{pT'(0) + (1 -p)r'(A)}, (9) 
X 2 = K{(l-p)r'(0)+pT'(-A)}, (10) 

A3 = ^{(i-p)r'(A)+pr'(-A)}, (11) 

with multiplicity l,Np—l,N(l—p) — l and 1 respectively. 
T'(x) is defined as -^T(x). Ao, which vanishes identically, 
always exists due to the invariance of Eq. (4) with respect 
to a constant shift of ipA an d "0s- ^i an( A ^2 are asso- 
ciated with the fluctuations of the individual oscillators 
belonging to the A-cluster and B-cluster, respectively. 
A3 corresponds to the fluctuation in A. Figure [| displays 
the eigenvalues versus A obtained using (g) with K = 1, 
which shows that all two-cluster states are unstable. It is 
important to note that (p, A') and (1 —p, A") are saddles 



which have negative eigenvalues of A2 and A 3 . (p, A'"), 
however, has a positive A3, which can be verified by the 
property A 3 cx j£p(A). 

Paradoxically, the system converges to unstable solu- 
tions. This counterintuitive fact may be understood if we 
assume that the pair of saddles (p, A') and (1 —p, A") are 
connected heteroclinically |i|,[l(| . All numerical results in 
Ref. [0 support this assumption. Although the hetero- 
clinicity is considered structurally unstable, this does not 
seem to apply to the particular class of systems under 
consideration. In the next section, it will be confirmed 
that (p, A') and (1 — p, A") are in fact connected hete- 
roclinically through an invariant subspaces, and argued 
how this structure is maintained stably. 



III. STRUCTURE OF THE HETEROCLINIC 
LOOP 

We first note a particular symmetry of our phase model 
given by Eq. (1) which is expressed as 



for all 



(12) 



V-i(i)=^(t) 



The above equation tells when the phases of some os- 
cillators are found identical at some time, they will obey 
completely the same dynamics thereafter or, equivalently, 
a point-cluster remains to be a point-cluster forever. This 
property will turn out crucial to the formation of hetero- 
clinic loops. 

We now argue how the phase model (|lj) can form hete- 
roclinic loops which is structurally stable. A generaliza- 
tion to the vector model of limit cycle oscillators will be 
given in the subsequent sections. Existence of a hetero- 
clinic loop connecting (p, A') and (1 — p, A") is clear if 
the following properties are satisfied: 

(X) (p, A') is a global attractor of W u {\ - p, A"), 

(Y) (1 - p, A") is a global attractor of W u (p, A'), 

where W u (p, A) represents the unstable manifold of 
(p, A). We work with an (N— l)-dimcnsional phase space 
throughout, by which the degree of freedom associated 
with a rigid rotation of the entire system is ignored. For 
an aid to the understanding of a little complicated sit- 
uation, it would be appropriate to display in advance a 
schematic picture of the heteroclinic loop under consid- 
eration in Fig. ||, where A-, A'/ and A'/' (i — 1,2,3) are 
the eigenvalues of (p, A'), (1 — p, A") and (p 7 A"') respec- 
tively. The definition of E x (X=A,B,AB) will be given 
later. 

Our argument will be based on the assumptions that 
there are three two-cluster states in the range p m i n < p < 
1 — Pmin and that the eigenvalues associated with these 
solutions satisfy certain stability properties. Specifically, 
the assumptions may be summarized as follows: 
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FIG. 5. Schematic representation of the structure of a 
heteroclinic loop, (p, A') and (1 —p, A") are the attractors of 
the invariant subspace Ea U -Bab and Eb U Eab respectively. 



(a) (p,A'), A") and (p,A'") exist, 

(b) A; > 0, 

(c) A!, < 0, A 3 < 0, 

(d) \'{ > 0, 

(e) Ag < 0, A 3 ' < 0, 

(f) A3" > 0, 

(g) Ai of all two-cluster states are positive, 

(h) ro = 0) > 0. 

Here, we assumed that both Np and N(l — p) are 
larger than 1 so that three independent eigenvalues A; 
(i = 1,2,3) exist. Note that the whole assumptions are 
satisfied by Eq. (|l|) with the coupling function given by 
Eq. (||) as stated in Sec. ||. 

Consider the tangent space around (p, A'), (p, A') has 
Np — 1 degenerate eigenvalues given by A^. Thus, the 
corresponding eigenvectors span a (Np — l)-dimensional 
unstable eigenspace which is denoted by Eb. Similarly, 
the eigenvectors corresponding to X' 2 span a (N(l — p) — 
l)-dimensional stable eigenspace which is denoted by 
Ea- An eigenvector corresponding to A3 spans the 1- 
dimensional stable eigenspace -Eab- In particular, the 
following statements hold. 

(b') Eb is the unstable subspace of (p, A'), 

X[ corresponds to the fluctuations which occur in the 
A-cluster (i.e. the phase-advanced cluster). Thus, the 
eigenspace Eb is associated with the disintegration of 
the A-cluster, while the B-cluster remains to be a point 
cluster. Similarly, the B-cluster is disintegrated in the 
eigenspace Ea, while the A-cluster remains to be a point 
cluster there. In the space -Eab, in contrast, these clus- 
ters remain to be point-clusters while their mutual dis- 
tance changes. Since a point-cluster must remain to be 
a point-cluster at any time, as noted at the beginning of 



this section, the space Eb U -Eab , on which the B-cluster 
is a point-cluster, gives an invariant subspace of dimen- 
sion Np. Similarly, Ea U Eab, on which the A-cluster 
is a point-cluster, is an invariant subspace of dimension 
N(l - p). Note that the unstable manifold W u (p, A') 
must coincide with Eb in the vicinity of (p, A'). This 
fact combined with the obvious fact that Eb is included 
in the invariant subspace Eb UEab leads to the following 
statements which hold globally: 

(xl) W u (p, A') is included by the invariant subspace 
E B UE AB . 

Arguments parallel to the above can be developed 
around (1 — p, A"), i.e., the state where the B-cluster 
is phase-advanced by A". From the assumed property 
(e), Eb UEab is the stable subspace of (1 —p, A"), which 
can be restated as follows: 

(e') (1 — p, A") is an attractor of E B U Eab- 

Therefore, if 

(x2) (1 — p, A") is a unique attractor of Eb U Eab, 
then we may assume 

(x3) (1 — p, A") is a global attractor of Eb U Eab- 

(xl) and (x3) give the sufficient conditions for (X) to 
hold. Similar discussion can be developed to derive (Y) 
via the assumption 

(y2) (p, A') is a unique attractor of Ea U Eab- 

Note that (p, A'") gives a source point lying between 
(p, A') and (1 —p, A") in the phase space as displayed in 
fig- I 

There still remains a problem about the validity of 
assumptions (x2) and (y2). It is sufficient to consider 
(x2) only. If the type of the attractors is limited to the 
two-cluster state, then it is obvious that (1 — p, A") is a 
unique attractor of Eb U Eab , as can be confirmed by the 
property (g). How about the possibility for a one-cluster 
state, namely, perfect synchrony, to become an attractor? 
The eigenvalues of the one-cluster state are N — 1-fold de- 
generate and given by KT'(x = 0), which is positive by 
the assumption (h) so that there is no stable manifold of 
one-cluster state. How about the stability structure of 
n(> 3)-cluster states? They could be attractors of the 
invariant subspaces with the same reason as for the two- 
cluster states even if they are unstable solutions. In the 
case of three or more clusters, however, the resulting het- 
eroclinicity would be even more complicated. Numerical 
simulations in the previous section, however, display the 
simple heteroclinicity between pairs of two-cluster solu- 
tions, implying the validity of the assumption (x2) and 
(y2) in Eq. (1) with Eq. (2). 

Convergence to the heteroclinic loop can be discussed 
similarly to the case of the heteroclinic orbit in a two- 
dimensional phase space, which was discussed in Ref. 
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JlO| . The result is that the system which is initially close 
to a heteroclinic loop converges to it provided 



rao 



< 1. 



(13) 



If this condition is satisfied, the heteroclinic loop is at- 
tracting. In numerical simulations, this convergence is 
established in a finite time due to the round-off errors. 

Although we have assumed the conditions (a)-(h) so 
far, our discussion is expected not to rely so heavily on 
the specific form of T(x). In fact, these conditions may be 
satisfied for a broader class of the coupling function. For 
instance, they are satisfied if we assume a simple shape 
of the coupling function such that r(x) decreases in the 
range — it < x < 0, while it increases otherwise (see Fig. 
^J). The reason is the following. The corresponding shape 
of p( A) turns out similar to the curve in fig. |[ so that we 
can define p m ; n similarly. For a given p satisfying p m - m < 
p < 1 — p m i ri , we obtain three states (p, A'), (p, —A") = 
(1 -p, A") and (p, A'"). Ai of each (p, A > 0) is positive 
because r'(0 < x < it) > 0. We can verify A3 < 0, 
A3 < and A 3 " > through the property that A3 is 
proportional to ^p(A). The one-cluster state turns out 
unstable since T'(0) > 0. Hence we have confirmed (a)- 
(h) except for A^A^' < 0. For the last properties to be 
satisfied, we need one more assumption, that is, r'(0) 
is not so large as to admit a region of p where both A2 
and X'2 are negative. Such region is expressed by p* < 
p < 1 — p*, where p* satisfies p mm < P* < 0.5. Then, 
via the assumptions (x2) and (y2), we obtain sufficient 
conditions for the existence of a heteroclinic loop between 
(p, A') and (1 — p, A") within the range p* < p < 1 - p*. 

In the previous works [^|JTo|l , the shape of the employed 
coupling functions satisfied this property. Our model dis- 
cussed in the next section also fulfills this condition in 
the weak coupling limit where the phase description is 
valid. Thus, we may regard the coupling function with 
this property of the shape as a typical class which admits 
heteroclinic loops. 

The existence of the phase space structure yielding a 
heteroclinic loop has thus been confirmed, which can be 
summarized as follows. For a given coupling function 
r(x), we can easily verify whether the conditions (a)- 
(h) are satisfied. Among these conditions, (a)-(f) con- 
stitute a necessary condition for the existence of a het- 
eroclinic loop between (p, A') and (1 — p, A"), while (g) 
and (h) support the assumptions (x2) and (y2). One 
may also consider the case where the roles of Ai and A2 
are reversed. The saddle connections are stably formed 
through the invariant subspace which exists for the sym- 
metry of equations of motion, or (|l2|). Thus, we conclude 
that the heteroclinic loop is robust under such small per- 
turbations that maintain the homogeneity of the popu- 
lation and the symmetry of the global coupling. 
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FIG. 6. Typical shape of the coupling function which 
admits a heteroclinic loop. This shape lead to conditions 
(a)-(h) except for X' 2 ,X'2 < 0. The last property also hold if 
r'(0) is not too large. 



IV. COUPLED LIMIT CYCLE OSCILLATORS 

Our argument on the existence and structural stability 
of the heteroclinic loop developed in the preceding sec- 
tion was based on the phase model with some assumed 
properties of the phase-coupling function. In this sec- 
tion, we discuss a specific coupled oscillator model in 
which heteroclinic loops are formed. From the model, 
a phase-coupling function of the desired properties for 
the existence of heteroclinic loops is derived through the 
method of the phase reduction. To our knowledge, the 
existence of the heteroclinic loop associated with slow 
switching has never been reported for vector models of 
oscillators. 

Consider a general system of coupled oscillators which 
are identical and all-to-all coupled: 



j t X z (t) = F(X l ) + ^J2 G ( X - X 



.))■ 



(14) 



Here Xi,F and G are m-dimensional real vectors and 
AT is a positive constant. Note that Eq. (14) satisfies the 
condition 



dt 



X ( (t)=X,(t) 



for all i, j, 



(15) 



which is similar to Eq. (12). Suppose that the local dy- 
namics is two-dimensional, i.e. X = (x,y), and the spe- 
cific forms of F and G are given by 



F x 



3xi 2 - x, ( 3 + y 4 - n 



1 — 5X; 



V, 



(16) 



GlXi, X j 



G, 



(17) 
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FIG. 7. Time series of the order parameter. The order 
parameter O is conveniently defined in the following way. Let 
tj (j = 0, 1, 2, ■ • •) denote the time at which the representative 
point of the j-th oscillator crosses a given section E in the 
1-oscillator phase space. The order parameter at time t — tjv 



EN 
i=1 exp 



is defined as 0(t = ijv) = 

generalization of Eq. (3). Since the oscillators cross E again 
and again, the order parameter at discrete times t = thN 
(k = 1, 2, • • ■) may be defined similarly. Note that 0(t) = 1 
when the oscillators are synchronized perfectly and 0(t) = 
when their phases are uniformly distributed. 

The corresponding equation X = F is called the 
Hindmarsh-Rose model [ pd{ which was originally pro- 
posed for a neural oscillator. Without coupling, i.e. 
K = 0, each unit becomes oscillatory if —11.5 < fi < 0.8 
fl3|| . The coupling is assumed to be diffusive, and in 
terms of neurophysiology, this corresponds to the electri- 
cal synapse formed by gap junctions |T^| . 

The parameter values are set to K = 0.1, N = 100 and 
fx = — 1. The intrinsic frequency then becomes u> ~ 1.0. 
We choose random initial conditions. Some numerical 
results obtained are summarized as follows. The system 
converges after a long transient to a two-cluster state 
which is periodic in time. Figure [?] displays a time se- 
ries of the order parameter. The relative population of 
the clusters generally depends on the initial condition. 
Convergence to the two-cluster state does not imply its 
stability, and is rather due to numerical artifacts. Actu- 
ally, when very small perturbations are given to the os- 
cillators independently, the clusters start to disintegrate, 
implying their linear instability. Such behavior is com- 
pletely similar to that of the phase model when hetero- 
clinic loops exist. We now show that the phase reduction 
of the above model produces a phase-coupling function 
which admits, based on the argument of the preceding 
section, heteroclinic loops. 

Coupled oscillators can be reduced to the phase model 
(Q) when the coupling is sufficiently weak M. There is 
a general formula for the phase-coupling function, and 
for a given dynamical-system model, this can be com- 
puted numerically. We did this for Eqs. (0), ( pif ) and 
(|l7|). The coupling function T(x) obtained is displayed in 



Fig. H which shows a typical shape admitting heteroclinic 
loops (see Fig. ||). Two-cluster solutions were sought and 
their stability analysis was done through Eqs. (j^)-(^TJ) . 
Then, we confirmed that the reduced model satisfies the 
conditions (a)-(h) for the existence of a heteroclinic loop 
and also the condition (O) for its stability. 



V. STRUCTURE OF THE HETEROCLINIC LOOP 
IN VECTOR MODELS 

The preceding arguments clarify the nature of the hete- 
roclinic loop in the framework of the phase model. In this 
section, we show that such arguments can be generalized 
to the original model for coupled limit-cycle oscillators in 
the vector form. 

It would be appropriate to start reconsidering the 
model given by Eqs. (|lj), © and @. Under suitable 
initial conditions, we obtain various two-cluster states. 
They correspond to the solutions (p, A) with A3 < in 
the phase-coupling limit, and these two-cluster states will 
be denoted by p-states. Also, the clusters corresponding 
to the phase-advanced and retarded clusters will respec- 
tively be called the A-cluster and B-cluster. For a given 
p-state, mN Lyapunov eigenvalues can be computed nu- 
merically, where m = 2 for the model under consider- 
ation. They can be classified into four groups Ao, A]S, 
A2S and A3S, and they will degenerate respectively into 
Xi (i — 0,1,2,3) in the phase-coupling limit. Note that 
each of Ais and A2S is composed of m eigenvalues, while 
A3S is composed of 2m — 1 eigenvalues. Ais corresponds 
to the fluctuations within the A-cluster, and each eigen- 
value of this group is Np — 1-fold degenerate. Similarly, 
A2S corresponds to the fluctuations within the B-cluster, 
and each eigenvalue there is N(l — p) — 1-fold degener- 
ate. A3S is associated with the relative motion between 
the clusters. Ao is identical to 0, resulting from the time- 
periodicity of the solutions. The maximum value of A^s 
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FIG. 8. Coupling function of the reduced model. The 
minimum of F(x) appears at a negative x. Such a shape of Y 
is typical when a heteroclinic loop exists. 
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FIG. 9. Lyapunov spectrum plotted against p. From the 
figure, it can be read out that p m in — 0.42. Solid lines show 
the eigenvalues obtained by the phase reduction, with excel- 
lent agreement with the spectrum of the original system. 

(z = 1,2,3) is denoted by Af ax , and their numerical val- 
ues are displayed in Fig. 0. 

The argument on the structure of a heteroclinic loop 
connecting the p-state and (1 — p)-state can be devel- 



III. We have 



oped quite in parallel with that in Sec. 
now to work with the (mN — l)-dimensional phase space, 
whereby we employ a surface of section to remove the ir- 
relevant degrees of freedom corresponding to the steady 
rotation of the whole population. The eigenspaces -Ea 
and E-q associated with Ais of p-state and (1 — p)-state, 
are now m{(N(l — p) — l}-dimensional and m(Np — 1)- 
dimensional, respectively, while the eigenspace £ab as- 
sociated with A3S is (2m — l)-dimcnsional. Then the 
argument in Sec. Ill still holds if we replace Aj with 
A™. Note that we assume the existence of an unstable 
state corresponding to (p, A'") which is hard to obtain 
numerically. 

The eigenvalues A™ ax are the ones which should co- 
incide with Xi in the phase-oscillator limit. As the cou- 
pling becomes stronger, the heteroclinic loop can persist 
as far as the existence- and stability properties of the 
two point-clusters are unchanged. Generally speaking, 
the stronger coupling makes point clusters more stable. 
In Fig. 9, this effect is already sizable for K = 0.1. As 
K becomes O(l), the two-cluster state gives way to a 
one-cluster state by which the heteroclinic loop disap- 
pears. We expect in general that the heteroclinic loop 
disappears when the coupling is so strong that the phase 
description completely breaks down. 



VI. APPEARANCE OF HETEROCLINIC LOOPS 
THROUGH DELAY-INDUCED BIFURCATIONS 

In globally coupled identical oscillators, one-cluster 
state is the easiest state to appear. This can be illus- 
trated by the following form of coupling: 



G(X i ,X j ) = X j (t)-X i (t). (18) 

Assuming Eqs. (^), ( [To] ) and (JlS|) , we obtain a stable 
one-cluster state even if K is very small. We may gener- 
ally expect that the heteroclinic loop cannot exist when 
the one-cluster state is stable. In the above model, it can 
be shown that the time-delayed coupling causes instabil- 
ity of the one-cluster state, which at the same time is 
accompanied by the appearance of the heteroclinic loop. 
The corresponding bifurcation is the so-called transcriti- 
cal bifurcation. 

Consider uniformly delayed coupling: 



G(Xi,Xj) = Xj(t - t) - Xi(i) 



(19) 



where r denotes the delay. Note that the symmetry prop- 
erty ( |l2| ) still holds when the coupling involves a uniform 
delay in the form of Eq. (|l9| ) . We will show some numer- 
ical results obtained for the system given by Eqs. (14), 
(16) and (19) where the parameter values are the same 
as in Sec. |v|. Without delay, the system under various 
initial conditions immediately converges to a one-cluster 
state. As r is increased, the one-cluster state persists up 
to a critical value r c beyond which the cluster splits into 
two and at the same time heteroclinic loops are formed. 
In this case of the parameters, this critical value is about 
0.18. 

This result can be understood by a phase reduction of 
our model which is applicable when the coupling is weak. 
The reduced model takes the form 



d K N 

-Vi(t)=w + -^r[^(t)-^-(*- 



% (20) 



where w ~ 1.0 at fi = — 1. Since the second term on the 
right-hand side is much smaller than the first term by 
assumption, (EG) is further reduced to the form 






x 



n 



FIG. 10. Coupling functions. The solid line is obtained 
from (|l|), ([l6|) and @ with r = 0, while the dotted line 
is obtained just by a phase shift of solid line by —ujt. The 
effect of the delay is equivalent to a simple modification of the 
coupling function in the weak-coupling limit. The modified 
coupling function is a typical shape admitting a heteroclinic 
loop provided the condition r > xq/uj is satisfied. 



7 



d K 

-A(t) = w + - T Mt) - + <t]. (21) 

i=i 

Thus, there is no explicit delay in coupling, while its effect 
has now been converted to a phase shift of the coupling 
function by lot. The situation is illustrated if Fig. 10. 
The stability of a one-cluster state depends entirely on 
the sign of T'(u>r). Thus, the one-cluster state is stable 
for small r which admits T'(u>t) < 0. As r is increased, 
the one-cluster state becomes less stable, and at t — xq /lj 
it becomes unstable where xq is defined as the value of x 
which minimizes T(x). Note that r c obtained numerically 
would agree with Xq/oj(~ 0.13) for sufficiently small K. 
For r > xq/lj, the coupling function assumes a typical 
shape under which a heteroclinic loop exists. At r = 0.3, 
for example, the conditions (a)-(h) and ( |l3| ) are satisfied 
in this reduced model. This transition occurs through a 
transcritical bifurcation. 

Let x denote the mutual phase difference between the 
clusters. The equation obeyed by x can be derived simi- 
larly to the derivation of Eq. (6), and takes the form 

j t x{t) = K{{2p - l)r(wr) + (1 - p)T(x + ujt) 

- P T(-x + wr)}. (22) 

The above equation has a trivial solution x = 0. For 
small x, the right-hand side can be expanded in powers 
of x. We find that provided p ^ 0.5 the right-hand side 
involves x 2 term as the lowest nonlinearity. This means 
that the trivial solution loses stability via a transcritical 
bifurcation. This occurs at r = xq/uj. We also find that 
as r is increased a saddle-node bifurcation occurs slightly 
before the transcritical bifurcation. For r > xq/ui, the 
stable manifolds associated with the saddle-node and 
transcritical bifurcations connect smoothly and forms a 



.... heteroclinic 




X<Tc Tc T<T C 



FIG. 11. Schematic representation of the bifurcation 
structure. The trivial solution x — loses its stability at 
t — t c via a transcritical bifurcation. Two solid lines existing 
for r > r c correspond to (p, A') and (1 — p, A") respectively, 
each being unstable with respect to the y- or z-direction al- 
ternately. A heteroclinic loop is formed between these two 
solutions as explained in Sec. |ril| 



heteroclinic loop. Let y (resp. z) denote a certain di- 
rection taken on E-q (resp. Ea)- How the bifurcation 
in question occurs is explained schematically in Fig. [ll]. 
The same structure of bifurcation leading to the forma- 
tion of heteroclinic loops holds for some range of p where 
X' 2 and X'2 are both negative. Note that at p = 0.5 the 
term of x 2 vanishes so that a pitchfork bifurcation occurs 
at r = xq/lu and no saddle- node bifurcation occurs be- 
fore. A heteroclinic loop is formed similarly to the case 
of p f 0.5. 

In coupled oscillators, the appearance of the hetero- 
clinic loop may seem to be pathological. The results in 
this section, however, imply that the heteroclinic loop 
appears in a broad class of weakly coupled oscillators 
including those studied so far provided the uniformly de- 
layed coupling is introduced. Assuming a simple coupling 
such as Eq. (18), we find that when the assemblies are 
composed of relaxation oscillators, their phase coupling 
function is often characterized by a curve which is sharply 
decreasing in a small region while it is gradually increas- 
ing otherwise. Thus, the heteroclinic loop is expected to 
arise in the homogeneous assemblies of relaxation oscil- 
lators. 



VII. SLOW SWITCHING 

When the system involves heteroclinic loops, it ex- 
hibits a remarkable dynamics when perturbed weakly. 
In the analysis of a model of the form of Eq. (1), Hansel 
et al. JTo| applied week noise independently to the indi- 
vidual oscillators and observed the appearance of a very 
long time scale depending on the noise intensity (see Fig. 
|l|). Since the time scale here is associated with the al- 
ternation between two collective states (i.e. a pair of 
two-cluster states), they called this characteristic behav- 
ior of the system slow switching, and gave a successful 
explanation for it in terms of a weakly perturbed hete- 
roclinic loop. Their explanation may be summarized as 
follows. If a heteroclinic loop is attracting, i.e. if Eq. ( fl3| ) 
is satisfied, then the system approaches one saddle and 
then to the other alternately. Without noise, the minimal 
distance from each saddle should decrease exponentially 
in time. With noise, however, this distance will fluctu- 
ate but remain finite typically within the order cr, the 
square-root of the variance of the noise. In any case, the 
system stays for the most time close to one saddle or the 
other, so that the dynamics could be characterized dom- 
inantly by the local properties around the saddles. The 
time-interval T for a stay near a saddle may be estimated 
as 

T~--!-lno-, (23) 
X u 

where X u is the eigenvalue of the most unstable direc- 
tion. Thus, the period of the switching is logarithmically 
dependent on the noise intensity. 



8 



0(f) 



0.5 




2000 



4000 



6000 



t 



FIG. 12. Time series of the order parameter for a noisy 
system. The solid line shows slow switching, where a new 
time scale of dynamics appears. 

By including noise, Eq. (14) is generalized as 
d N 

-X^t) = F(Xi) +KY J G{X %1 Xj) + (24) 



noise, a slight violation of the symmetry property jl5| ) 
is expected to cause the same effect. Imagine a partic- 
ular case where a hetcroclinic loop is present under the 
symmetry condition (|l5|). The system is now perturbed 
slightly so that the condition ( |l5| ) is slightly violated. 
Assume that a pair of saddles, between which a hetero- 
clinic loop is formed in the symmetric case, still exists. 
Although a genuine heteroclinic loop could no longer ex- 
ist in the asymmetric system, the unstable manifold of 
one saddle will come close to the other saddle. If they 
are sufficiently close in the phase space, the situation is 
quite similar to the case of applied weak noise, leading 
to slow-switching. 

As an example, let fi in Eq. ( |l6| ) be Gaussian- 
distributed with variance p 2 . Numerical simulations ac- 
tually show slow switching without noise. Figure dis- 
playing the period versus /?, shows again a logarithmic 
law. The steepness of |T/lnp| is estimated to be 20, 
which implies that the gap e is now the order of p. Sim- 
ilar results are obtained when the other parameters, e.g. 
the delay r, are randomly distributed. We can also break 
the uniformity in the coupling, and consider a slightly 
random oscillator network. 



where £j is Gaussian white noise with variance 1 and the 
parameter a, assumed to be small below, indicates the 
intensity of noise. Let us now consider specific forms for 
F and G given by Eq. (16) and Eq. (19), respectively, 
where the parameter values are the same as before. Some 
results obtained numerically are the following. For < 
t < r c , the oscillators localize in a small range of phase, 
which we call noisy one-cluster state. Figure [l^ displays 
a time trace of the order parameter O, where r = 0.1 
and p = 10~ 7 . It is seen that O stays near 1. For r > 
r c , this highly coherent cluster becomes unstable, and 
O starts to oscillate. After a long transient, the system 
comes to exhibit slow switching between a pair of two- 
cluster states. For the most time the system stays close 
to one of the noisy two-cluster state, which is followed by 
a short period of cluster disintegration, then converging 
to another two-cluster state. This is demonstrated in Fig. 
(14) for r = 0.3. It is clear that the collective dynamics 
is then characterized by a new time scale corresponding 
to this slow switching. We define the period of switching 
T as the average time between the two successive local 
minima of O sufficiently after the transient. Logarithmic 
dependence of T on the noise intensity is clear from Fig. 
|l3| . The steepness of the T versus In cr curve after linear 
fitting is estimated as \T / ln<j| ~ 20, which suggests that 
A™ ax of this delayed coupling model is about 0.05. A™ ax 
can be easily estimated by Eq. (9) with the phase-shifted 
coupling function displayed in Fig. 10. We obtain Ai ~ 
0.065 at p = 0.5, which is close to the above estimation 
of A™ ax ~ 0.05. Note that the amplitude effect makes 
A™ smaller than Ai similarly to the case of Fig. 9. 

Slow-switching is thus the fate of the system when the 
heteroclinic loop is perturbed weakly. Besides external 



VIII. CONCLUSIONS 

In the present paper, the physical relevance of the spe- 
cific model employed has little been discussed. We used 
the Hindmarsh-Rose oscillator model which was origi- 
nally proposed as a model for neuronal oscillators. We 
have seen that the heteroclinic loop appears rather easily 
in an assembly of relaxation oscillators, so that something 
similar may be expected for homogeneous neuronal as- 
semblies subject to a constant external current. However, 
we did not intend modeling realistic neuronal populations 
with the Hindmarsh-Rose model, but it was used just for 
its simplicity in demonstrating the formation of the het- 
eroclinic loop in coupled oscillators especially when the 
coupling involves the time-delay. It could be understood 




FIG. 13. Switching period versus noise intensity. The line 
is a linear fitting of the data. 
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FIG. 14. Switching period versus the square-root of the 
variance of the parameter. The line is a linear fitting of the 
data. 



strong perturbation and its frequency becomes compa- 
rable with the intrinsic frequency of the oscillators, the 
dynamics would become even more complicated due to 
the nonlinear coupling between these modes of motion of 
comparable time-scales. Statistical mechanical approach 
to this problem would be interesting. As far as we have 
analyzed, however, slow switching vanishes well before its 
frequency comes close to the intrinsic frequency. It would 
be also interesting to find out coupled-oscillator models 
which have more robust structure of slow switching with 
respect to the symmetry-breaking perturbations. 
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from our whole argument that the same conclusion would 
be valid for a much broader class of coupled oscillator 
models. In particular, any population model of coupled 
oscillators would be capable to form the heteroclinic loop 
provided that it satisfies the symmetry condition (15). 
If we associate the switching phenomenon with biolog- 
ical rhythms, the novel properties of this phenomenon 
would yield lots of suggestive ideas. One of the intrigu- 
ing subjects there is rhythm splitting In connection 
with neural networks, a more realistic model should be 
considered, and a study in this direction is now under 
progress. 

It may appear that the heteroclinic loop in question 
is something which could not go beyond some mathe- 
matical curiosity because the symmetry property ( |l5| ) on 
which it crucially relies would be more or less violated 
in real systems. However, the associated phenomenon 
of slow switching seems to be of much greater physical 
relevance, because the strict symmetry need not be re- 
quired there. Since the noise, heterogeneity and delay 
are commonplace in macroscopic systems, some indica- 
tion of slow switching could well be detected in the real 
world. In the case of mechanical oscillators, for instance, 
it would not be difficult to obtain oscillators which are 
almost identical. Global coupling might also be realized 
through electric circuit jL5| , vibrating board |l^] , surface 
motion of water |lj] and so on. A certain class of sur- 
face chemical reactions under oscillatory conditions may 
provide globally coupled identical oscillators. 

The characteristic frequency of slow switching will be- 
come shorter when the symmetry breaking perturbations 
become stronger. At the same time, the amplitude of the 
oscillating order parameter become smaller as the pertur- 
bations become stronger. The switching phenomenon is 
expected to vanish when the strength of the perturbation 
exceeds a critical value after which the order parameter 
of the system becomes stationary. Realistic examples 
of slow switching, if any, would correspond to the case 
of strong perturbations. If slow switching survives the 



* Email: kori@ton.scphys.kyoto-u.ac.jp 
[1] S. Watanabe and S.H. Strogatz, Physica D 74, 197 
(1994). 

[2] K. Wiesenfeld, P. Colet, and S.H. Strogatz, Phys. Rev. 

Lett. 76, 404 (1996). 
[3] D. Hansel, G. Mato, and C. Meunierm, Europhys. Lett. 

23, 367 (1993). 
[4] K. Okuda, Physica D 63, 424 (1993). 
[5] Y. Kuramoto, Chemical Oscillation, Waves, and Turbu- 
lence (Springer, New York, 1984). 
[6] R. R. Aliev and V. N. Biktashev, J. Phys. Chem. 98, 

9676 (1994). 
[7] A. T. Winfree, J. Theor. Biol. 16, 15 (1967). 
[8] J. Buck, Q. Rev. Biol. 63, 265 (1988). 
[9] KR. Delaney, A. Gelperin, M. S. Fee, J. A. Flores, 

R. Gervais, D. W. Tank, and D. Kleinfeld, Proc. Natl. 

Acad. Sci. USA 91, 669 (1994). 
[10] D. Hansel, G. Mato, and C. Meunier, Phys. Rev. E 48, 

3470 (1993). 

[11] J. L. Hindmarsh and R. M. Rose, Proc. R. Soc. London 
B 221, 87 (1984). 

[12] L. Abbott and E. Marder, in Method in Neuronal Mod- 
eling, edited by C. Koch and I. Segev (MIT Press, Cam- 
bridge, 1989). 

[13] S. Raghavachari and J. A. Glazier, Phys. Rev. Lett. 82, 
2991 (1999). 

[14] H. Iglesia, J. Meyer, A. Carpino and W. Schwartz, Sci- 
ence 290, 799 (2000). 

[15] K. Wiesenfeld, P. Colet, and S.H. Strogatz, Phys. Rev. 
Lett. 76, 404 (1996). 

[16] S.H. Strogatz and I. Stewart, Sci. Am. 269, No. 6, 102 
(1993). 

[17] K. Yoshikawa, N. Oyama, M. Syoji, and S. Nakata, Am. 
J. Phys. 59, 137 (1991). 



10 



